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Abstract 

Context. The instrument BLAST (Balloon-borne Large-Aperture Submillimeter Telescope) performed the first deep and wide extra- 
galactic survey at 250, 350 and 500 /jm. The extragalactic number counts at these wavelengths are important constraints for modeling 
the evolution of infrared galaxies. 

Aims. We estimate the extragalactic number counts in the BLAST data, which allow a comparison with the results of the P(D) analysis 
of Patanchon et al. (2009). 

Methods. We use three methods to identify the submillimeter sources. 1) Blind extraction using an algorithm when the observed field 
is confusion-limited and another one when the observed field is instrumental-noise-limited. The photometry is computed with a new 
simple and quick point spread function (PSF) fitting routine (FASTPHOT). We use Monte-Carlo simulations (addition of artificial 
sources) to characterize the efficiency of this extraction, and correct the flux boosting and the Eddington bias. 2) Extraction using a 
prior. We use the Spitzer 24 pm galaxies as a prior to probe slightly fainter submillimeter flux densities. 3) A stacking analysis of the 
Spitzer 24 pm galaxies in the BLAST data to probe the peak of the differential submillimeter counts. 

Results. With the blind extraction, we reach 97 mjy, 83 mjy and 76 mjy at 250 pm, 350 pm and 500 pm respectively with a 95% 
completeness. With the prior extraction, we reach 76 mJy, 63 mJy, 49 mjy at 250 pm, 350 pm and 500 pm respectively. With the 
stacking analysis, we reach 6.2 mJy, 5.2 mJy and 3.5 mjy at 250 pm, 350 pm and 500 pm respectively. The differential submillimeter 
number counts are derived, and start showing a turnover at flux densities decreasing with increasing wavelength. 
Conclusions. There is a very good agreement with the P(D) analysis of Patanchon et al. (2009). At bright fluxes (>100 mjy), the 
Lagache et al. (2004) and Le Borgne et al. (2009) models slightly overestimate the observed counts, but the data agree very well 
near the peak of the differential number counts. Models predict that the galaxy populations probed at the peak are likely z ~ 1.8 
ultra-luminous infrared galaxies. 

Key words. Cosmology: observations - Galaxies: statistics - Galaxies: evolution - Galaxies: photometry - Infrared: galaxies 



1. Introduction 

Galaxy number counts, a measurement of the source surface 
density as a function of flux density, are used to evaluate the 
global evolutionary photometric properties of a population 
observed at a given wavelength. These photometric properties 
mainly depend on the source redshift distribution, spectral 
energy distribution (SED), and luminosity distribution in a 
degenerate way for a given wavelength. Even though this 
is a rather simple tool, measurements of number counts at 
different observed wavelengths greatly help in constraining 
those degenera c ies. B a ckward evoluti o n mod e ls, among these 
Charv & Elba zl d2001l); ILagache et all d2004l) : iGrappioni et al, 
d2005|) ; iFranceschini et al] d2009t): iLe Borgne et al. 
i200%: iPearson & Khanl d2009t) : iRowan-Robinsonl d2009): 



ValianteetalJ ([2009) are able to broadly reproduce (with 
different degrees of accuracy) the observed number counts from 
the near-infrared to the millimeter spectral ranges, in addition 
to other current constraints, like such as measured luminosity 
functions and the spectral en ergy distribution o f the Cosmic 
Infrare d Background (CI B) 

lHauser et al] Il998t ILagache et 



their redshift distributions, or the m ean temperature or col ors of 
galaxies, as is shown for instance in lLeFloc'hetaD(l2009l) from 
Spitzer 24 pm deep observations. 

One key spectral range lacks valuable data to get accurate 
constraints as yet: the sub-millimeter range, between 160 pm 
and 850 pm, where some surveys were conducted on small 
areas. Fortunately this spectral domai n is intensively stu died 
with the BLAST balloon experiment dDevlin etal]|2009l) and 
the Herschel and Planck space telescopes. This range, al- 
though it is beyond the maximum of the CIB's SED in wave- 
length, allows us to constrain the poorly-known cold compo- 
nent of galaxy SED at a redshift greater than a few tenths. 
Pioneering works h ave measured the local luminosity function 
dDunne et al]|2000T) a nd shown that most milli-Jansky sources 
lie at redshifts z > 2 dlvison et alJl2002b IChapman et alJ l2003a. 
120051: llvison et all 120051: iPope et al]|2005l 120061) . Other works 
showe d that the galaxies SED selected in the submillimeter 



1996; 
1999; 



Fixsen et al 



Gispert et al 



(Puget et al. 

HIT 

2000; lHauser & Dwekl 120011: iKashlinskvl 120051: ILagache et al 
2005; iDole et al] 12006). In the details, however, the models 
disagree in some aspects like the relative evolution of luminous 
and ultra-luminous infrared galaxies (LIRG and ULIRG) and 



range dBenford et al.lll999t IChapman et al.lE003bT: ISaiina et all 
2003 1: iLewis et al .1120051: iBeelen et alJl2006t . iKovacs et al]|2006: 
Saiin a et al] 12006; Michalow ski et al] 120091) can have typically 
warmer temperatures and higher luminosities than galaxies se- 
lected at other infrared wavelengths. 

Data in the submillimeter wavelength with increased sen- 
sitivity are thus needed to match the depth of infrared sur- 
veys, conducted by S pitzer in the mid - and far-infrared with 
the MIPS instrument dRieke et al]l2004l) at 24 pm, 70 pm and 
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160 pm (ICharv et al.ll2004t iMarleau et al.ll2004t iPaoovich et al 



2004; Dole et al.ll2004l 



Fraver et al. 2006a b; Rodiahiero et al 
2006; I Shupe et al.1 120081; iFraver et al.1 l2009t iLeFloc'h etal 
2009; IBethermin et al.l l2010h as well as the near-infrared 



range with the IR AC instrument dFazio et al. 2004b) between 
3.6 pm and 8.0 um dFazio et al.ll20 04a; Franceschini et al. 2006; 
Sulli van et all l2007t iBarmbv et alJ 120081: iMagdis et al J 12008; 
Ashbv et all 120091) . Infrared surveys have allowed the resolu- 
tion of the CIB by identifying the contributing sources - directly 
at 24 and 70 um, or indirectly tro ugh stacking at 160 pm 
dDole et al.ll2006l IBethermin et alj|2010h . 

Although large surveys cannot solve by themselves all the 
unknowns about the submillimeter SED of galaxies, the con- 
straints given by the number counts can greatly help in unveiling 
the statistical SED shape of submillimeter galaxies as well as the 
origin of the submillimeter background. 

The instrument BL AST (Balloon - borne Large-Aperture 
Submillimeter Telescope, Pascal e et al.l (12008)) perf ormed the 



first w ide and deep survey in the 250-500 pm range (Devlin et al 
2009) before the forthcoming Herschel results. Marsd en et al 



(2009) show that sources detected by Spitzer at 24 um emit the 
main part of the submillimeter background. iKhan et a I] (120091) 
claimed that only 20% of the CIB is resolved by the sources 
brighter than 17 mJy at 350 ^m. lPatanchon et al] d2009l) has per- 
formed a P{D) fluctuation analysis to determine the counts at 
BLAST wavelength (250 pm, 350 pm and 500 pm). In this pa- 
per we propose another method to estimate the number counts 
at these wavelengths and compare the results with those of 
iPatanchon etall d2009h . 



2. Data 

2.1. BLAST sub-millimeter public data in the Chandra Deep 
Field South (CDFS) 

The BLAST holds a bolometer array, which is the precursor of 
the spectral and photometric imaging receiver (SPIRE) instru- 
ment on Herschel, at the focus of a 1.8 m diameter telescope. It 
observes at 250 fim, 35 um and 500 um , with a 36", 42" and 
60" beam, respectively (Truc h et alj|2009l) . 

An observation of the Chandra Deep Field South (CDFS) 
was performed during a long duration flight in Antarctica in 
2006, and the data of the two surveys are now p ublic: a 8.7 deg 2 
shallow field and a 0.7 deg 2 confusion-limited dDole et alJ 2004) 
field in the center part of the first one. We use the non-beam- 
smoothed maps and associated point spread function (PSF) dis- 
tributed on the BLAST websit^B The signal and noise maps 
were generated by the SANEPIC algorithm (Pata nchon et al] 
2008)." 



2.2. Spitzer 24 um data in the CDFS 

Several infrared observations were performed in the CDFS. 
The Spitzer Wide-Field InfraRed Extragalactic (SWIRE) sur- 
vey overlaps the CDFS BLAST field at wavelengths between 
3.6 u m an d 160 pm. We used only the 24 pm band, which is 
80% complete at 250/Jy. The completeness is defined as the 
probability to find a source of a given flux in a catalog. The Far- 
Infrared Deep Extragalactic Legacy (FIDEL) survey is deeper 
but narrower (about 0.25 deg 2 ) than SWIR E and 80% com plete 
at 57 pjy at 24 pm. We used the Betherm in et al] (Hp 10) cata- 
logs constructed from these two surveys. These catalogs were 




Figure 1. Position of sources brighter than the 95% complete- 
ness flux at 250 pm in deep zone. Top: initial map. Bottom: 
residual map. The area out of the mask are represented darker. 
These l°xl° map are centered on the coordinates (RA,Dec) = 
(3h32min30s,-27°50')- The horizontal axis is aligned with the 
right ascension. 



http://www.blastexperiment.info 



extracted with SExtractor dBertin & Arnoutsl 19961) and the pho- 
tometry w as performed with the allstar routine of the DAOPHOT 
package (Stetson 1987). The completeness of this catalog was 
characterized with Monte-Carlo simulations (artificial sources 
added on the initial map and extracted). 

3. Blind source extraction and number counts 

We started with a blind source extraction in the BLAST bands. 
Each wavelength was treated separately. For each wavelength 
we defined two masks: a shallow zone (about 8.2 deg 2 ) covering 
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Figure 2. Extragalactic number counts at 250 fim in the BLAST 
data, diamond: counts deduced from the source catalog on the 
whole shallow field; square: counts deduced from the catalog of 
the deep part of the field; triangle: counts deduced from catalog 
of the deep part of the field with a 24 /jm prior (this measurement 
gives only a lower limits to the counts); cross: counts computed 
with a stacking analysis; grey asterisk: counts computed with a 
P(D) analysis b y Patanc hon et alj (120091) ; grey short dashed line: 
Laga che et al.l (12004]) model predictio n; grey long dashed line 
and grey area: iLe B orgne et al. (2009) model prediction and 1- 
o~ confidence area. 



the whole field except the noisier edge; and a deep zone (about 
0.45 deg 2 ) in the center of the confusion-limited area. We used 
different extraction methods in the shallow zone and the deep 
one, but the photometry and the corrections of the extraction bias 
were the same. 



3.1. Detector noise-limited extraction (shallow zone) 

In the shallow zone we used the non-smoothed map and the cor- 
responding map of the standard deviation of the noise. The map 
was then cross-correlated by the PSF. The result of this cross- 
correlation is 



vO'o, jo) 



+N 

z 

i=-N j=-N 



m(i + i, jo + j) x PSF(i, j), 



(1) 



where m com ,(io, jo) is the flux density in the pixel (;'o, jo) of the 
cross-correlated map, m(i, j) the flux density in the pixel (z, j) 
of the map, and PSF(i, j) the value of the normalized PSF in 
the pixel (i, j) (the center of the PSF is in the center of the pixel 
(0, 0)). The PSF size is (2N+l)x(2N+l) pixels. The standard de- 
viation of the noise in the cross-correlated map is thus 



riconviio, jo) 



|] nHio + iJo+j)xPSFmj), 



E 

i=-N j=-N 



(2) 



where n (n conv ) is the initial (cross-correlated) map of the stan- 
dard deviation of the noise. 

We found the pixels where m conv /n com , > 3 and kept the local 
maxima. The precise center of the detected sources was com- 
puted by a centroid algorithm. This low threshold caused lots 
of spurious detections, but helped to deblend the fluxes of 3 to 
4-sigma sources and avoided to overestimate their fluxes. We 
could thus limit the flux boosting effect. A final cut in flux af- 
ter the PSF fitting photometry eliminated the main part of these 
sources. We performed the extraction algorithm on the flipped 
map (initial map multiplied by a factor of -1) to check it. We 
found few spurious sources brighter than the final cut in flux de- 
termined in the Sect. 13.41 We found a spurious rate of 12%, 11% 
and 25% at 250 fim, 350 yum and 500 yum, respectively. 

3.2. Confusion-limited extraction (deep zone) 

In the confusion-limited zone we also used a non-smoothed map. 
In this region the noise is dominated by the confusion and not by 
the instrumental noise. Consequently, the method based on in- 
strumental noise presented in the Sect. 13. II is not r e levant. We 
used an atrou wavelet filtering (Starck et ail 1 19991 iDole et all 
|200T1) to remove fluctuations at scales larger than 150". Then we 
divided the resulting map by o-fn, ere( i map , which is the standard 
deviation of the pixel values on the filtered map in the working 
area. We finally kept local maxima with a signal greater than 3. 
The center of the sources was also determined by a centroid al- 
gorithm. The initial map and the cleaned map are shown in Fig. 
Q] When we flip the map, we find no spurious source brighter 
than the final cut in flux determined in Sect. 13.41 



3.3. A simple and quick PSF fitting routine: FASTPHOT 

For both noise- and confusion-limited extraction, we apply the 
same quick and simple PSF fitting routine on the non-beam- 
smoothed map. This routine fits all the detected sources at the 
same time and is consequently efficient for deblending (although 
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no source was blended in this case; but source-blending will be 
an issue for an extraction using a prior, detailed in Sect. [4] We 
suppose that the noise is Gaussian and the position of sources is 
known. We then maximize the likelihood 



L(m\S) = Y\ C(n)xexp 



pixels 



2n 2 



(3) 



where m and n are the map and the noise map. PS F Xhy . is a unit- 
flux PSF centered at the position (x,,j i ), which are the coordi- 
nates of the i-th source. These coordinates are not necessarily 
integers. C(n) is a normalization constant and depends only of 
the value of the noise map. S is a vector containing the flux of 
the sources. 

The value of S, which maximizes the likelihood, satisfies the 
following linear equation stating that the derivative of the likeli- 
hood logarithm equals zero 



Vi, : 



dlog(L(m\S)) 



dSi 



AS +B, 



where A is a matrix and B a vector defined by 

/'.S /\ ,,. x/'.V/\ , v 

n 2 

pixels 



A = K) = - 2 



B = (bt) = 2 

pixels 



PS F x , y . x map 



(4) 



(5) 



(6) 



To perform this operation fast, we used a 70"x70" (respectively 
90"x90" and 110"xll0") PSF at 250 yum (respectively 350 fxm 
and 500 fim). This PSF, provided by the BLAST team, is the 
response for a unit-flux source and takes into account all the fil- 
tering effects. We used the conjugate gradient method to solve 
this equation quickly. 

This routine was tested with 200x200 pixels simulated maps 
containing 400 sources at a known positions with a beam of 10 
pixels FWHM. The flux of all sources was perfectly recovered in 
the case where no noise was added. This routine (FASTPHOT) 
performs simultaneous PSF fitting photometry of 1000 sources 
in less than 1 second. It is publicly availably 

3.4. Completeness and photometric accuracy 

The completeness is the probability to detect a source of a given 
flux density. We measured it with a Monte-Carlo simulation. We 
added artificial point sources (based on PSF) on the initial map at 
random positions and performed the same source extraction and 
photometry algorithm as for the real data. A source was consid- 
ered to be detected if there was a detection in a 20" radius around 
the center of the source. TableQ]gives the 95% completeness flux 
density (for which 95% of sources at this flux are detected) for 
different wavelengths and depths. 

The photometric noise was estimated with the scatter of the 
recovered fluxes of artificial sources. We computed the standard 
deviation of the difference between input and output flux. This 
measurement includes instrumental and confusion noise (cr t0 , = 

yfrfnstr + ^conf)- ^ ne resu l ts m ' e given in Table Q] In the deep 
area, the photometric uncertainties are thus dominated by the 
confusion noise. The estimations of the confusion noise between 



on the IAS website http://www.ias.u-psud.fr/irgalaxies/ 



the deep and shallow areas are consistent. It shows the accuracy 
and the consistency of our method. 

No te that the uncertainties on flux densities in the lDve et al.l 
(2009) catalog (based only on instrumental noise) are conse- 
quently largely underestimated in the confusion-limited area. 
Indeed, their 5 <x detection threshold (based only on instrumen- 
tal noise) at 500 /im in the deep zone corresponds to 1 .76 cr if 
we also include the confusion noise. 

The faint flux densities are overestimated due to the classical 
flux boosting effect, . This bias was measured for all bands for 
60 flux densities between 10 mJy and 3 Jy with the results of the 
Monte-Carlo simulations. The measured fluxes were deboosted 
with this relation. We cut the catalogs at the 95% completeness 
flux, where the boosting factor is at the order of 10%. Below 
this cut, the boosting effect increases too quickly to be safely 
corrected. We also observed a little underestimation at high flux 
of 1%, 0.5% and 0.5% at 250 fjm, 350 fim and 500 yum. It is 
due to FASTPHOT, which assumes that the position is perfectly 
known, which is not true, especially for a blind extraction. 

3.5. Number counts 

We computed number counts with catalogs corrected for boost- 
ing. For each flux density bin we subtracted the number of spuri- 
ous detections estimated in the Sects. [3~Tl and l3.2l to the number 
of detected sources and divided the number of sources by the 
size of the bin, the size of the field and the completeness. 

We also applied a corrective factor for the Eddington bias. 
We assumed a distribution of flux densities in dN/dS oc S~ r 
with r — 3 ± 0.5. This range of possib l e valu es for r was es- 
timated considering the [Patanchon et al. (2009) counts and the 
iLaeache et alJ (2004) and Le Borg ne et alJ ([2009) model predic- 
tions. We then randomly kept sources with a probability given 
by the completeness and added a random Gaussian noise to sim- 
ulate photometric noise. Finally we computed the ratio between 
the input and output number of sources in each bin. We applied 
a correction computed for r — 3 to each point. We estimated the 
uncertainty on this correction with the difference between cor- 
rections computed for r — 2.5 and r = 3.5. This uncertainty was 
quadratically combined with a Poissonian uncertainty (cluster- 
ing effects are negligible due to the little number of sources in 
the map, see appendix lAl. 

The calibration uncertainty of BLAST is 10%, 12% and 13% 
at 250 fim, 350 fim and 500 /mi respectively dTruch et alj |2009). 
This uncertainty is combined with other uncertainties on the 
counts. The results are plotted in Fig. [2] and given in Table [2] 
and interpreted in Sect. [6] 

3.6. Validation 

We used simulations to valid ate our method. We g enerated 50 
mock catalogs based on the Patancho n et al ] (120091) counts, and 
which covered 1 deg 2 each. These sources are spatially homo- 
geneously distributed. We then generated the associated maps at 
250 //m. We used the instrumental PSF, and added a gaussian 
noise with the same standard deviation as in the deepest part of 
real map. 

We performed an extraction of sources and computed the 
number counts with the method used in the confusion lim- 
ited part of the field (Sect. 13 ,2b . We then compared the output 
counts with the initial counts (Fig. [3]). We used two flux density 
bins: 100-141 mJy and 141-200 mJy. We found no significant 
bias. The correlation between the two bins is 0.46. The neigh- 
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20.3 
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131 76 
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26.4 


17.6 


16.7 


16.5 



Table 1. 95% completeness flux density and photometric noise for different depths at different wavelengths. The instrumental 
noise is given by the noise map. The total photometric noise includes the instrumental and confusion noise and is determined by 

Monte-Carlo simulations. The confusion noise is computed with the formula cr co „f = Jo~^ ot - o~\ nstr - 



bor po ints are thus not anti-correlated as in the iPatanchon et alJ 
(2009) P(D) analysis. 

The same ve ri ficatio n was done on 20 
Fernandez-Conde e t"al] d2008f) simulations (based on the 



Lagache et alJ (2004) model). These simulations include cluster- 
ing. This model overestimates the number of the bright sources, 
and the confusion noise is thus stronger. The 95% completeness 
is then reach at 200 mJy. But there is also a very good agreement 
between input and output counts in bins brighter than 200 mJy. 
We found a correlation between two first bins of 0.27. 
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Figure 3. Number counts at 250 yum deduced from a blind 
extraction for 50 real izations of a simulation based on the 
Pata nchon et al.l ((2009) counts. The horizontal solid line repre- 
sents the input count value. The lower panel is the result of the 
100-141 mJy bin, the upper panel is the 141-200 mJy bin. 



4. Source extraction using Spitzer 24 pm catalog as 
a prior 

In addition to blind source extraction in the BLAST data 
(Sect. [3} we also performed a source extraction using a prior. 

4.1. PSF fitting photometry at the position of the Spitzer 
24 pm 

The catalogs of infrared galaxies detected by Spitzer contain 
more sources than the BLAST catalog. The 24 /im Spitzer PSF 
has a Full Width at Half Maximum (FWHM) of 6.6". It is smaller 
than the BLAST PSF (36" at 250 /an). Consequently, the po- 



sition of the Spitzer sources is known with sufficient accuracy 
when correlating with the BLAST data. 

We applied the FASTPHOT routin e (Sect. l33Tl at the posi- 
tions of 24 fim sources. We used the Bethermin et al. I J20T0I) 
SWIRE catalog cut at S 24 = 250 pJy (80% completeness). In 
order to avoid software instabilities, we kept in our analysis only 
the brightest Spitzer source in a 20" radius area (corresponding 
to 2 BLAST pixels). The corresponding surface density is 0.38, 
0.49 and 0.89 Spitzer source per beamj at 250 pm, 350 pm and 
500 pm, respectively. 

This method works only if there is no astrometrical offset be- 
tween the input 24 pm catalog and the BLAST map. We stacked 
the BLAST sub-map centered on the brightest sources of the 
SWIRE catalog and measured the centroid of the resulting arti- 
ficial source. We found an offset of less than 1". It is negligible 
compared to the PSF FWHM (36" at 250 pm). 

We worked only in the central region of the deep confusion- 
limited field (same mask as for blind extraction), where the pho- 
tometric noise is low. 



4.2. Relevance of using Spitzer 24 pm catalog as a prior 

The 5 25o/5 24 (S 350/5 24 or S 500 /S 24) c °l° r ls not constant, and 
some sources with a high color ratio could have been missed in 
the prior catalog ( especi ally high-redshift starbursts ). We used 
the lLagacheeF al. (2004) and lLeBorgne et all d2009l) models to 
estimate the fraction of sources missed. We selected the sources 
in the sub-mm flux density bin and computed t he 24 pm flux 
densi ty distribution (see Fig. According to the Lagach e et al.l 
(2004) model, 99.6%, 96.4% and 96.9% of the sub-mm selected 
sources 4 are brighter than S 24 = 250 yuJy f or a selection at 
250 pm , 350 pm and 500 pm, respectively The lLe Borgne et al] 
(120091) model gives 99.8%, 98.3% and 95.0%, respectively. 



4.3. Photometric accuracy 

The photometric accuracy was estimated with Monte-Carlo arti- 
ficial sources. We added five sources of a given flux at random 
positions on the original map and add them to the 24 pm cata- 
log. We then performed a PSF fitting and compared the input and 
output flux. We did this 100 times per tested flux for 10 flux den- 
sities (between 10 and 100 mJy). In this simulation we assumed 
that the position of the sources is exactly known. It is a reason- 
able hypothesis due to the 24 yum PSF FWHM (6.6") compared 
to the BLAST one (36" at 250 yum). 

We did not detect any boosting effect for faint flux densities 
as expected in this case of detection using a prior. For a blind 
extraction there is a bias of selection toward sources located on 
peaks of (instrumental or confusion) noise. This is not the case 



3 the beam solid angles are taken as 0.39 arcmin 2 , 0.50 arcmin 2 and 
0.92 arcmin 2 at 250 pm, 350 pm and 500 yum respectively. 
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4.5. Sub-mm/24 color 
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Figure 4. Flux density distributi on at 24 pm of the 54 mJy < 
S 350 < 83 mJy sour ces. The Lagache et al. (2004) model is plot- 
ted in black and the Le Borg ne et alj ( 2009 ) model is plotted in 



grey. The dashed line represents the cut of our catalog at 24 pm. 

for an extraction using a prior, for which the selection is per- 
formed at another wavelength. 

The scatter of output flux densities is the same for all the in- 
put flux densities. We found a photometric noise cr s of 2 1 .5 mJy, 
18.3 mJy and 16.6 mJy at 250 pm, 350 pm and 500 /mi, respec- 
tively. It is slightly lower than for the blind extraction, for which 
the position of source is not initially known. 

4.4. Estimation of the number counts 



From the catalog described in Sect. 14.11 we give an estimation 
of the submillimeter number counts at flux densities fainter than 
reached by the blind-extracted catalog. We cut the prior catalog 
at 3 os , corresponding to 64 mJy, 54 mJy and 49 mJy at 250 /im, 
350 pm and 500 pm, respectively. We worked in a single flux 
density bin, which is defined to be between this value and the cut 
of the blind-extracted catalog There is no flux boosting effect, 
but we needed to correct the Eddington bias. The completeness 
could not be defined in the same way as for the blind extrac- 
tion, because the selection was performed at another wavelength. 
We thus cannot suppose power-law counts, because the selection 
function is unknown and the distribution of the extracted sources 
cannot be computed. 

The Eddington bias was estimated with another method. We 
took the sub-mm flux of each of the sources selected at 24 /mi 
and computed how many sources lie in our count bin. We added 
a Gaussian noise <x s to the flux of each source to simulate the 
photometric errors. We computed the number of sources in the 
counts bin for the new fluxes. We then compute the mean of the 
ratio between the input and output number of sources in the se- 
lected bin for 1000 realizations. The estimated ratios are 0.42, 
0.35 and 0.21 at 250 /mi, 350 /mi and 500 /mi, respectively. 
These low values indicate that on average the photometric noise 
introduces an excess of faint sources in our flux bin. This ef- 
fect is strong because of the steep slope of the number counts, 
implying more fainter sources than brighter sources. The results 
are interpreted in the Sect. [6] 



4 the bins are defined as 64 to 97 at 250 fim, 54 to 83 at 350 fim and 
49 to 76 at 500 fim 



In this part we work only on S > 5 cr s sources of the cata- 
log described in Sect. 14.11 to avoid bias due to the Eddington 
bias in our selection. At 250 fim, we have two sources verify- 
ing this criterion with a S250/S24 color of 16 and 60. No sources 
are brighter than 5 c ry at larger wav elengths. For this cut in flux 
(S250 > 5 cry), the lLagache et al.l d2004l) and iLe Borgne et al.l 
(2009) models predict a mean S250/S24 color of 39 and 41, re- 
spectively. The two models predict a mean redshift of 0.8 for this 
selection, and the K-correction effect explains these high colors. 



5. Non-resolved source counts by stacking analysis 

5.1. Method 

In order to probe the non-resolved source counts, we used same 
method as Bethermin et al. (2010), i.e. the stacking analysis ap- 
plied to number counts (hereafter "stacking counts"). We first 
measured the mean flux at 250 fim, 350 fim or 500 fim as a func- 
tion of the 24 fim flux (5250, 350 or 500 = /OS 24))- This measure- 
ment was performed by stacking in several S 24 bins. We used the 
Bethermi n et alj d2010t) catalog at 24 fim of the FIDEL survey. 
It is deeper than the SWIRE one used in Sect. [4] but covers a 
smaller area (0.25 deg 2 ). The photometry of stacked images was 
performed with the PSF fitting method (Sect. l3.3l l, and the uncer- 
tainties on the m ean flux are computed with a bootstrap method 
dBavo uzet 2008). We then computed the counts in the sub-mm 
domain with the following formula: 



dN 



dS 



sub nun 



dN 



S submtn-f{S 24) 



dS 



24 



idS 



subrnm 



S24 



dS 



24 



(7) 



We show in appendix 151 that the clustering effect can be ne- 
glected. The results are given in Tableland are plotted in Fig. 

m 



5.2. Validity of the stacking analysis in the sub-mm range 

There are 1.8, 2.4 and 4.5 S24 >70/dy sources per BLAST beam 
at 250 fim, 350 fim and 5 00 um, respectively. We thus stacked 
several sources per beam. Bethermin et al. (2010) showed that 
the stacking analysis is valid at 160 fim in the Spitzer data, where 
the size of the beam is similar to the BLAST one. 

To test the validity of the stacking analysis in the BLAST 
data from a Spitzer 24 /mi catalog, we generated a simula- 
tion of a 0.25 deg 2 with a Gaussian noise at the same level 
as for the real map and w ith source clustering, following 
iFernandez-Conde et alJ d2008l) . We stacked the 24 pm simulated 
sources per flux bin in the BLAST simulated maps. We measured 
the mean BLAST flux for each 24 pm bin with the same method 
as applied on the real data. At the same time we computed the 
mean sub-mm flux for the same selection from the mock cata- 
log associated to the simulation. We finally compared the mean 
BLAST fluxes measured by stacking with those directly derived 
from the mock catalog to estimate the possible biases (see Fig. 
[5). The stacking measurements and expected values agree within 
the error bars. We notice a weak trend of overestimation of the 
stacked fluxes at low flux density (S24 < 200 pJy) however, but 
it is still within the error bars. We can thus stack 24 pm Spitzer 
sources in the BLAST map. 
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Table 2. Number counts deduced from source extraction. The not normalized counts can be obtained dividing the S 25 .dN/dS 
column by the S mean column. 
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high 24 yum flux densities (S 24 > 400 yuJy) the derivative is al- 
most constant and small (<20 and compatible with zero), mean- 
ing that the observed sub-mm flux density does not vary much 
with S 24. For these flux bins we select only local sources and do 
not expect a strong evolution of the color. At fainter 24 fim flux 
densities the observed decrease can be explained by redshift and 
K-correction effects, as above. 

The color in the faintest 24 fim flux density bin (70 to 102 
fiJy) is slightly fainter than in the neighboring points. It can be 
due to the slight incompleteness of the 24 fim catalog (about 
15%), which varies spatially across the field: the sources close 
to the brightest sources at 24 yum are hardly extracted. The con- 
sequence is a bias to the lower surface density regions, leading 
to a slight underestimation of the stacked flux measurement. 



Figure 5. Ratio between the mean flux density at 250 fim found 
by the stacking anal ysis and the expected flux for d ifferent S24 
bins. It is based on a iFernandez-Conde et alj ((2008) simulation 
of a 0.25 deg 2 field with a noise and PSF similar to the BLAST 
deep region. 

5.3. Mean 24 fim to sub-mm color deduced by stacking 
analysis 

The stacking analysis allowed to measure the mean 24 yum to 
sub-mm colors of undetected sub-mm galaxies. These colors 
depends on the SED of galaxies (or K-correction) and the red- 
shift distribution in a degenerate way. The S su i, mm /S 24 color and 
dS submm/dS 24 as a function of S 24 are plotted in Fig. [6] 

The colors are higher for the fainter 24 fim flux (S 24 < 
100 yuJy). This behavior agrees with the model expectations: the 
faint sources at 24 fim lie at a higher mean redshift than the 
brighter ones. Due to the K-correction, the high-redshift sources 
have a brighter sub-mm/24 color than local ones. 

The colors found by the stacking analysis are lower than 
those obtained by an extraction at 250 yum (Sect. 14.5b . It is an 
effect of selection. The mid-infrared is less affected by the K- 
correction than the sub-mm , and a selection at this wavelength 
selects lower redshift objects. We thus see lower colors because 
of the position of the SED peak (around 100 fim rest-frame). 

We also investigated the evolution of the derivative 
dS su bmm/dS 24 as a function of S 24, which explicits how the ob- 
served sub-mm flux increases with the 24 fim flux densities. At 



5.4. Accuracy of the stacking counts method on BLAST with 
a Spitzer 24 fim prior 

Beth ermin et al.l d2010h showed that the stacking counts could 
be biased: the color of sources can vary a lot as a function of the 
redshift. The assumption of a single color for a given S24 is not 
totally realistic and explains some biases. We used two simulated 
catalogs (containing for each source S 24, S 250, S350 and S 500) 
to est imate this effect: a first one based on the lLagache et al.l 
(2004) model that cove r ed 20 x2.9 deg 2 and a second one based 
on the iLe Borgne et all d2009) model and that covered 10 deg 2 . 
The large size of these simulations allows us to neglect cosmic 
variance. 

In order to compute the stacking counts, we first computed 
the counts at 24 fim from the mock catalog. Then we computed 
the mean S 250. 350 or 500 flux density (directly in the catalog) in 
several S 24 bins to simulate a stacking. We finally applied the 
Eq.Qto compute stacking counts at the BLAST wavelengths. 

The ratio between the stacking counts and the initial counts 
is plotted in Fig.[8]for the two mock catalogs. Between 1 mJy and 
10 mJy we observe an oscillating bias. This bias is less than 30% 
at 250 fim and 50% at other wavelengths. When the flux becomes 
brighter than 25 mJy at 250 fim (18 mJy at 350 fim and 7.5 mJy 
at 350 yum), we begin to strongly underestimate the counts. The 
analysis of real data also shows a very strong decrease in the 
counts around the same fluxes (see Fig . [TJ . Consequently, we cut 
our stacking analysis at these fluxes and we applied an additional 
uncertainty to the stacking counts of 30% at 250 fim (50% at 
350 fim and 500 fim). 
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Figure 6. Black solid line: dS submmldS 24 as a function of S 24. 
Grey dashed line: S su bmm/S 24 color as a function of S 24. 

Using the 24 /mi observations as a prior to stack in the 
BLAST bands seems to give less accurate results than in the 
Spitzer MIPS bands. For a given S 24 flux, the sub-mm emis- 
sion can vary a lot as a function of the redshift. But the sim- 
ulations shows that this method works for faint flux densities. 
It is due to the redshift selection which is similar for faint flux 
densities (see Fig. [9]) and very different at higher flux densities 
(see Fig.flOb. For example, S24 ~100 /dy sources are distributed 
around z — 1 .5 with a broad dispersion in redshift. S350 ~ 4 mJy 
(based on averaged colors, 4 mJy at 350 /mi corresponds to 
S24 ~ 100 /dy) sources have quite a similar redshift distribution 
except an excess for z > 2.6. At higher flux densities (around 
2 mJy at 24 /mi) the distribution is very different. The majority 
of the 24 //m-selected sources lies at z < 1 and the distribution 
of 350 /mi-selected sources peaks at z ~1.5. Another possible 
explanation is that fainter sources lies near z=l and are thus se- 
lected at the 12 /mi rest-frame, which is a very go od estimator of 
the inf rared bolometric luminosity according to Spinoglio et al.l 
(119951) . 

In order to limit the scatter of the sub-mm/24 color, 
we tried to cut our sample into two redshift boxes fol- 
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Table 3. Number counts deduced from stacking. The not nor- 
malized counts can be obtained dividing the S 2 5 .dN/dS column 
by the S column. 



lowing the iDevlin et ail d2009t) IRAC color criterion ([3.6]- 
[4.5]=0.068([5.8]-[8.0])-0.075). But we had not enough signal 
in the stacked images to perform the analysis. 



6. Interpretation 

6.1. Contribution to the CIB 

We integrated our counts assuming power-law behavior between 
our points. Our points are not independent (especially the stack- 
ing counts), and we thus combined errors linearly. The con- 
tribution of the individually detected sources (S250 > 64 mJy, 
S350 > 54 mJy, S500 > 49 mJy) is then 0.24+° |* nW.m^sr 1 , 
0.06!°'^ nW.m 2 .sr' and 0.01+°°| nW.m 2 .sr f at 250 /mi, 
350 fim and 500 ^m, re spectively. Considering the total CIB 
level of lFixsen et al.l (1 19981) (FIRAS absolute measurement), we 
resolved directly only 2.3%, 1.1% and 0.4% at 250 /mi, 350 /urn 
and 500 //m, respectively. 

The populations probed by the stacking counts 
(S250 > 6.2 mJy, S350 > 5.2 mJy, Sgoo > 3.5 mJy) emit 
5.0!?^ nW.m 2 .sr', 2.8^ nW.mlsr 1 and 1.4+ 2 ' nW.m 2 .sr' 
at 250 /mi, 350 /mi and 500 //m, respectively. This corresponds 
to about 50% of the CIB at these three wavelengths. 

6.2. Comparison with \Patanchon et al\ Ij200$) 

The agreement between our r esolved counts built fr om the cata- 
logs and the P(D) analysis of Pata nchon et al.l (120091) is excellent 
(Fig. |2]i. We confirm the efficiency of the P(D) analysis to recover 
number counts without extracting sources. The stacking counts 
probe the flux densities between 6 mJy and 25 mJy at 250 //m 
(between 5 mJy and 13 mJy at 350 /mi and 3 mJy and 7 mJy 
at 500 fim). In this range there is only one P(D) point. At the 
three BLAST wavelengths the P(D) points agree with our stack- 
ing counts (Fig . [2]). Our results thus confirm the measurement of 
iPatanchon et alJ d2009l) and give a better sampling in flux. 
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Figure 7. Number counts at BLAST wave lengths coming from the data (points) and the models (lines). Short dashed line: initial 
(black) and stacking (grey) counts fr om the lLagache et al.l (2004) mock catalog; long dashed line: initial (black) and stacking (grey) 
counts from lLe Borgne et alJ (l2009t) : diamond: stacking counts built with the FIDEL catalog; square: stacking counts built with the 
SWIRE catalog; grey vertical dot-dash line: flux cut for stacking counts. 
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Figure 8. Ratio between stacking counts and initial counts for two mock catalogs; short dashed line: Lagache et al. ( 2004) catalog; 
long dashed line: iLe Borgne et alJ d2009t) ; grey vertical dot-dash line: flux cut for stacking counts; horizontal dot line: estimation of 
the uncertainty intrinsic to the stacking method. 



6.3. Comparison with ground-based observations 

We compared our results with sub-m m ground-based obser- 
vations of SHARC. iKhan et al.l d2007l) estimated a density of 
S350 > 13 mJy sources of 0.84+lg' arcmin~ 2 . For the same cut, 
we found 0.26+0.13 arcmin -1 , which agrees with their work. 
Our measurement (175± 75 sources. de g~ 2 brighter than 25 mJy) 
agrees also with that of ICoppinet al. (2008) ones at the same 
wavelength (200-500 sources. deg -2 brighter than 25 mJy). 

We also compared our results at 500 /an with the SCUBA 
ones at 450 fxm. iBorvs etall d2003l) find 140+^° gal.deg" 2 for 
S450 >100 mJy. We found 1.2+1.0 gal.deg -2 . We significantly 



disagree with them. IBorvs et alJ {2003) claim 5 4-cr detections 
in a 0.046 deg 2 field in the Hubble deep field north (HDFN). 
These five sources are brighter than 100 mJy. We find no source 
brighter than 100 mJy in a 0.45 deg 2 field at 350 /mi nor at 
500 /urn. The cosmic variance alone thus cannot explain this 
difference. A possible explanation is that they underestimated 
the noise level and their detections are dominated by spurious 
sources. It could also be due to a calibration shift (by more than 
a factor 2). The observation of the HDFN by Herschel will al- 
low us to determine whether that these bright sources might be 
spurious detections. 
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Figure 9. Solid line: Distribution in redshift of the sources with 
102 yuJy < S 24 < 150 /Jy for the mock ca talogs generated 
with th e ILagache et all d2004 (black) and the|Le Borgne etall 
(2009?) (grey) models; Dashed line: Distribution in redshift of 
the sources with 4 mJy < S 250 < 6 mJy (determined using the 
mean 250/24 color). 
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Figure 10. Solid line: Distribution in redshift of the sources 
with 2213 yJ y < 5 74 < 3249 /Jy for the moc k catalogs gener- 
ated w ith the ILagache et al.l ((2004) (black) and lLe Borgne et all 
(2009j) (grey) models; Dashed line: Distribution in redshift of 
sources with 47 mJy < S 250 < 62 mJy (determined using the 
mean 250/24 color). 



We also compared our results with the estimations based 
on lensed sources a t 450 /mi with SCUBA dSmail et alj|2002t 
iKnudsen et~ai1l2006l) . For example, iKnudsen et all (120061) find 
2000-50000 sources. deg 2 brighter than 6 mJy. It agrees with our 
3500+™° sources.deg 2 . 

6.4. Comparison with the\Laaache etali $2004) and 
\Le Boraneei al. (2009) models 

At 250 //m and 350 //m the measured resolved source 
coun ts are signifi c antly lowe r (by about a factor of 2) than 
the ILagache et all d2004 and iLe Borgne et all d2009l) models. 



Figure 11. Mean redshift ( solid line) of source s for differ- 
ent fluxes at 350 fo r the ILagache et al.l (120041) (black) and 
ILe Borgne et al.l (12009") (grey) models and corresponding in- 
frared luminosity defined in Sect. \6.5\ ( dashed line). 



Nevertheless, our cou nts are within the confidence area of 
Le Borgn e et alJ (120091) . The same effect ( models overestimat- 
ing the counts) was observed at 160 fjm dFraver et al.l 12009: 
Bethermin et al It indicates that the galaxies' SED or the 

luminosity functions used in both models might have to be re- 
visited. At 500 /mi our counts and both models agree very well, 
but our uncertainties are large, which renders any discrimination 
difficult. 

Concerning the stacking counts, they agree very well with 
the two models. Nevertheless, our uncertainties are larger than 
30%. We thus cannot check if the disagreement observed be- 
tween the Lagache et al.1 (120041) m odel and the stacking counts 
at 160 /im (iBethermin et alJ2010h of 30% at S m = 20 mJy still 
holds at 250 /mi. 

6.5. Implications for the probed populations and the models 

We showed that the two models nicely reproduce the sub-mm 
counts, especially below 100 mJy. We can thus use them to esti- 
mate which populations are constrained by our counts. For each 
flux density bin we computed the mean redshift of the selected 
galaxies in both models. We then used the SEDs given by the 
models at that mean redshift and at that flux bin and derived the 
infrared bolometric luminosity. The luminosities are shown in 
dashed lines in Fig.QT]for 350 //m as an example, and the red- 
shift is given in solid lines. 

The stacking counts reach 6.2 mJy, 5.3 mJy and 3.5 mJy at 
250 fim, 350 /mi and 500 /mi, respectively. This corresponds to 
faint ULIRGs (L IR * 1.5 x 1O 12 L ) around z = 1.5, 1.8 and 2.1 
at 250 /mi, 350 fim and 500 /zm, respectively. Our measurements 
show that the predicted cold-dust emissions (between 100 /mi 
and 200 /mi rest frame) of this population in the models are be- 
lievable. 

At 250 //m and 350 /mi the resolved sources (S 350 > 
85 mJy) are essentially z ~ 1 ULIRGs (L m > 10 12 L o ) 
and HyLIRGs (L/r > 10 13 L©) according to the models. In 
Laga che et alJ d2004l) the local cold-dust sources contribute at 
ver y bright flux (> 200mJ y). This population is not present in 
the lLe Borgne et afl d2009l) model. It explains the difference be- 
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tween the two models for flu xes brighter than 100 mJy at 350 
fim (see Fig.Qj}- At 500 ^m. lLagache et al.1 (12004 predict that 
bright counts are domi nated by local cold-dust populations and 
iLe Borgne et alJ ([2009) that they are dominated by medium red- 
shift HyLIRGs. Nevertheless, there is a disagreement with the 
observations for this flux density range, suggesting that there 
could be less HyLIRGs than predicted. But these models do not 
currently include any AGN contrib ution, which is small except at 
luminosities higher than 1 12 L Q dLacv et al.ll2004t iDaddi et alJ 
l2007HValiante et alJl2009h . 

7. Conclusion 

Our analysis p rovides new stacking co unts, which can be com- 
pared with the [Patanchon et al.l ([2009) P(D) analysis. We have 
a good agreement between the different methods. Nevertheless, 
some methods are more efficient in a given flux range. 

The blind extraction and the extraction using a prior give a 
better sampling in flux and slightly smaller error bars. The P(D) 
analysis uses only the pixel histogram and thus looses the in- 
formation on the shape of the sources. The blind extraction is a 
very efficient method for extracting the sources, but lots of cor- 
rections must be applied carefully. When the confusion noise to- 
tally dominates the instrumental noise, the former must be deter- 
mined accurately , and the catalog flux limit must take this noise 
dDole et al.ll2003l) into account. 

Estimating the counts from a catalog built using a prior is 
a good way to deal with the flux boosting effect. This method 
is based on assumptions however. We assume that all sources 
brighter than the flux cut at the studied wavelength are present 
in the catalog extracted using a prior. We also assume a flux dis- 
tribution at the studied wavelength for a selection at the prior 
wavelength to correct for the Eddington bias. Consequently an 
extraction using a prior must be used in a flux range where the 
blind extraction is too affected by the flux boosting to be accu- 
rately corrected. 

P(D) analysis and stacking counts estimate the counts at flux 
densities below the detection limit. These methods have differ- 
ent advantages. The P(D) analysis fits all the fluxes at the same 
time, where the stacking analysis flux depth depends on the prior 
catalog's depth (24 fim Spitzer for example). But the P(D) anal- 
ysis with a broken power-law model is dependent on the number 
and the positions of the flux nodes. Th e uncertainty due to th e 
parameterization was not evaluated by [Patanchon et al. (2009). 
The stacking counts on the other hand are affected by biases 
due to the color dispersion of the sources. The more the prior 
and stacked wavelength are correlated, the less biased are the 
counts. A way to overcome this bias would be to use a selection 
of sources (in redshift slices for example), which would reduce 
the color dispersion, and the induced bias; we did not use this 
approach here because of a low signal-to-noise ratio. 

The stacking and P(D) analysis are both affected by the clus- 
tering in different ways. For the stacking analysis this effect de- 
pends on the size of the PSF. This effect is small for BLAST 
and will be smaller for SPIRE. The clustering broadens the pixel 
histogram. Patanchon et al ] (120091) show that it is negligible for 
BLAST. Clustering will probably be an issue for SPIRE. The 
cirrus can also affect t he P(D) analysis and broaden the peak. 
Patancho n et al.l d2009l) use a high-pass filtering that reduces the 
influence of these large scale structures. 

The methods used in this paper will probably be useful to 
perform the analysis of the Herschel SPIRE data. The very high 
sensitivity and the large area covered will reduce the uncer- 
tainties and increase the depth of the resolved source counts. 



Never theless, according to the models (e.g. ILe Borgne et al.l 
(2009)), the data will also be quickly confusion-limited and it 
will be very hard to directly probe the break of the counts. The 
P(D) analysis of the deepest SPIRE fields will allow us to con- 
strain a model with more flux nodes and to better sample the 
peak of the normalized differential number counts. The instru- 
mental and confusion noise will be lower, and a stacking analysis 
per redshift slice will probably be possible. These analyses will 
give stringent constraints on the model of galaxies and finally on 
the evolution of the infrared galaxies. 
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The expected results for mean stacking of an N non-clustered 
populations is 



M(8) = S s xPSF(0) + 



JT'S- 



(B.l) 



where M is the map resulting from stacking, 8 the distance to 
the center of the cutout image, S , the mean flux of the stacked 
population. The integral is an approximation because the cen- 
tral source is treated in the first term. This approximation is to- 
tally justified in a strongly confused field where the number of 
sources is enormous. PSF is the instrumental response and is 
supposed to be invariant per rotation (8 = corresponds to the 
center of this PSF). ^ is the number of the source per flux unit 
and per pixel. We assume an absolute calibration. The integral in 
the equation lB.ll is equal to the CIB brightness 



(B.2) 



This term is constant for all pixels of the image and corresponds 
to a homogeneous background. 

The stacked sources can actually be autocorrelated. The 
probability density to find a stacked source in a given pixel and 
another in a second pixel separated by an angle 8 (p(0)) is linked 
with the angular autocorrelation function ((d(8)) by 



p(ff)=p I s (l+<o(8)), 



(B.3) 



where p s is the number density of the stacked source. 

If we assume that there is no correlation with other popu- 
lations, the results of the stacking of N autocorrelated sources 
is 



Appendix A: Effect of clustering on the 
uncertainties of number counts 



Bethermi n et al.l (1201 Oh showed how the clustering is linked 
wi th the uncertain t ies of the counts. We used the formalism 
of Bether min et al.l (1201 Ol) to estimate the effect of the cluster- 
ing on our BLAST counts. There are few sources detected at 
250 pm and the BLAST coverage is inhomogeneous. It is con- 
sequently very hard to estimate the clustering of the resolved 
population. We thus us ed the clustering measured at 160 pm by 
Bethermin et al. and assumed a 250/160 color equal to 

unity. We then used the same method to compute the uncertain- 
ties. We then compare the uncertainties with and without clus- 
tering. Neglecting the clustering implies an underestimation of 
the uncertainties on the counts of 35% in the 203-336 mJy bin 
at 250 pm, and less than 20% in the other bins. We can thus 
suppose a Poissonian behavior, knowing that the Poisson ap- 
proximation underestimates the error bars for the 203-336 mJy 
bin at 250 pm. Nevertheless, our model of clustering at 250 pm 
has strong assumptions (single 250/160 color, same clustering at 
250 pm as measured at 160 pm), and it would be more conser- 
vative to update it with Herschel clustering measurements. 



Appendix B: Effect of clustering on stacking 

B.1. A formalism to link clustering and stacking 

The clustering can bias the results o f a stacking. We present a 
formalism based on lBavouze i (l2008l) work. 



M{0) = S , x PSF(0) + / C/B , s (l + (0(6)) * PSF(0) + I CIB , ns , (B.4) 

where Icib,s and Icis.ns is the CIB contribution of stacked and 
non-stacked sources. If we subtract the constant background of 
the image, we find 



M(8) = S s x PSF(8) + I CIB , S x u){8) * PSF(8). 



(B.5) 



The second term of this equation corresponds to an excess of 
flux due to clustering. This signal is stronger in the center of 
the stacked image. The central source appears thus brighter than 
expected, because of the contribution due to clustering. 

The flux of the central stacked source computed by PSF- 
fitting photometry is 



JfMxPSFdD. _ 

S mes — "7. ~ ~~ — S s + S clusi 

j j PSF 2 dQ. 



(B.6) 



where S c i us , the overestimation of flux due to clustering is given 
by 



* clus 



IciB,s X 



/ / ((w * PS F) x PS F)d£l 
J J PSF 2 dQ. 



(B.7) 



Basically, the stronger the clustering, the larger the bias. In ad- 
dition, the wider the PSF, the larger the overestimation. The 
stacked signal can be dominated by the clustering, if the angu- 
lar resolution of the instrument is low compared to the surface 
density of the sources (like Planck, c.f. Fernand ez-Conde et al.l 
d2010h ) or if strongly clustered populations are stacked. 
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B.2. Estimation of the bias due to clustering 

The estimation of S c i„ s with Eq. lB.7l requires particular hypothe- 
ses. The stacked population is S24 > 70 /iJy sources detected 
by Spitzer. Their contribution to the CIB is 5.8 nW.m~ 2 .sr _1 , 
3.4 nW.m^.sr -1 and 1.4 nW.irT 2 .sr~' at 250 /mi, 350 //m and 
500 /mi, respectively (estimated by direct stacking of all the 
so urces). Following the clustering of 24 /mi sources estimated 
bv lBethermin et alTd2010l) . we suppose the following autocorre- 
lation function: 



The excess of flux due to clustering (S c i us ) is then 0.44 mJy, 
0.35 mJy and 0. 16 mJy at 250 /mi, 350 //m and 500 /mi, respec- 
tively. This is significantly lower than the bootstrap uncertainties 
on these fluxes. We can thus neglect the clustering. 

B.3. Measurement of the angular correlation function by 
stacking 

This new formalism provides a simple tool to measure the angu- 
lar autocorrelation function (ACF) from a source catalog. This 
method uses a map called "density map". One pixel of this map 
contains the number of sources centered on it. It is equivalent of 
a map of unit flux sources with the PS F — 6 (Dirac distribution). 
The result of the stacking is thus 




(B.8) 



M(0) = p s x 6(6) + p. s (l + oj(0j). 



(B.9) 



The ACF can then be easily computed with 



V6» + 0, oj{0) = 



M(6) 



- 1. 



(B.10) 



